Problem Set 4: Non-linear Models and GAMs

PLS 206 — Applied Multivariate Modeling

Author

Grey Monroe

Published

January 1, 2026

Note

Due: Monday, October 26, 1:10 PM (start of class) — submitted through Canvas.

Questions 1–2 use Monday’s lecture (curves from new columns, the quadratic and its vertex). Questions 3–10 use Wednesday’s (reading a GAM, choosing k, by =, and nls()). The week’s code page has working code for every step; the hints point to its parts and sections.

Submission Instructions

Upload two files to Canvas:

  1. Written answers — a .pdf document with your written responses, interpretations, and figures.
  2. R script — a .R file with all the code you used to produce your answers.

List the first and last names of any collaborators in both files.

Naming convention: HW4emailID.R (code) and HW4emailID.pdf (written answers), where emailID is the part of your UC Davis email before the @.

Example: Grey Monroe would submit HW4gmonroe.R and HW4gmonroe.pdf.

Feel free to use Google, StackOverflow, and other resources if you are stuck with code, and post questions on Slack if anything is unclear.

Every model you plot: show the raw data, the fitted curve, and its 95% confidence band, with the response on the y-axis.

The data

Three datasets, all read straight from the course website:

  • climate.csv: global temperature anomalies (Temp), one row per year from 1880 to 2005, plus columns for greenhouse gases, land use, ozone, aerosols, volcanic activity, Earth’s orbit and solar activity (Greenhouse_gas, Land_use, Ozone, Aerosols, Volcanic, Orbit, Solar), factors that could explain the temperature.
  • hemispheres.csv: temperature anomalies (Temp) for the Northern and Southern Hemispheres (Hemisphere), one row per hemisphere per year from 1880 to 2020.
  • leaf_light.csv: a simulated light-response experiment. 120 leaves (60 sun leaves and 60 shade leaves, leaf_type), each measured once at one light level (PAR, µmol photons m⁻² s⁻¹) for net photosynthesis (A, µmol CO₂ m⁻² s⁻¹).
library(ggplot2)
library(mgcv)
climate <- read.csv("https://greymonroe.github.io/PLS_206/data/climate.csv")
hemi <- read.csv("https://greymonroe.github.io/PLS_206/data/hemispheres.csv")
leaf <- read.csv("https://greymonroe.github.io/PLS_206/data/leaf_light.csv")

Part 1: Is the Warming a Straight Line?

Q1. Plot Temp against Year, with the fitted line from lm(Temp ~ Year) and its 95% confidence band. Then plot the model’s residuals against Year. Is a straight line adequate? Describe the pattern in the residuals.

resid(model) gives the residuals; add a horizontal line at zero with geom_hline(yintercept = 0). For the band, use the code page’s fit_ci() or predict(model, newdata, interval = "confidence"). Script Part A, sections 2 and 3; slide Residuals: always look.

Q2. Fit a quadratic, lm(Temp ~ poly(Year, 2, raw = TRUE)), and plot it with its 95% band. Calculate the vertex from the coefficients: in what year does the quadratic say temperature was lowest? Is the vertex a good estimate of the coolest period in these data? One or two sentences.

The vertex is -b1 / (2 * b2), with b <- coef(model), b1 = b[2], b2 = b[3]. To judge it, look at the coldest years in the data (climate[order(climate$Temp), ][1:5, ]) and at the shape of the data around the vertex: a parabola is symmetric. Script Part A, sections 5 and 6; slides Top of the hill and A parabola is symmetric.

Part 2: A GAM for the Warming

Q3. Fit gam(Temp ~ s(Year), data = climate, method = "REML"). From summary(), report the edf, the p-value of the smooth, and the deviance explained. In one sentence, say what the edf means here. What does the intercept estimate equal, and why?

Compare the intercept with mean(climate$Temp). Script Part D, section 19; slides The summary, line by line, The intercept is the mean of y and Smooths are centred at zero.

Q4. Run gam.check() on your Q3 model and report the k-index and its p-value. Refit with k = 20 and plot both fitted curves (with their bands) on the data. Did the curve change? Did the edf get close to its ceiling? Using these two answers, explain why the low k-index here is not a sign that k is too small, and suggest what it does point to.

gam(Temp ~ s(Year, k = 20), data = climate, method = "REML"). The k-index p-value comes from a random shuffle, so it changes slightly each run (set.seed(206) before gam.check() makes it repeatable). Think about the rule: double k, refit, and see whether the curve moves. Then ask whether one year’s residual tells you anything about the next year’s: cor(resid(m)[-1], resid(m)[-nrow(climate)]). Script Part D, section 21; slides The k check: double it and see and k = 8 and k = 16: the same curve.

Q5. Fit one GAM per predictor for Greenhouse_gas, Land_use and Solar (Temp ~ s(Greenhouse_gas), and so on) and compare their deviance explained. Which one tracks temperature best? Why can’t this comparison on its own show what causes the warming?

Deviance explained is on the last-but-one line of each summary(). For the second part, plot each predictor against Year. Script Part D, section 19.

Part 3: Is the Warming Even Across the Planet?

Q6. Make Hemisphere a factor, plot Temp against Year coloured by hemisphere, and fit gam(Temp ~ Hemisphere + s(Year, by = Hemisphere), data = hemi, method = "REML"). Why does the model need Hemisphere on its own, as well as inside by =? Use the HemisphereSouthern row of the summary in your answer.

hemi$Hemisphere <- factor(hemi$Hemisphere). Each smooth is centred on zero, so what is left to carry each hemisphere’s average level? Compare the two parametric coefficients with tapply(hemi$Temp, hemi$Hemisphere, mean). Script Part E, section 24; slide by = needs the factor too.

Q7. Plot both fitted curves with their 95% bands on the data. Use predict() to get each hemisphere’s fitted temperature in 1980 and in 2020, with 95% confidence intervals. Which hemisphere warmed more over those 40 years, and do the two hemispheres’ intervals in 2020 overlap?

nd <- expand.grid(Year = c(1980, 2020), Hemisphere = levels(hemi$Hemisphere)), then predict(m, nd, se.fit = TRUE); the interval is fit ± 1.96 * se.fit. For the plot, do the same on Year = 1880:2020. Script Part A, section 3, and Part D, section 20; slide For a figure: predict().

Part 4: Light-Response Curves with nls()

A light-response curve has parameters with a biological meaning: \(A = A_{max}\,(1 - e^{-\alpha\,\mathrm{PAR} / A_{max}}) - R_d\), where \(A_{max}\) is the maximum gross photosynthesis (so net photosynthesis levels off at \(A_{max} - R_d\)), \(\alpha\) is the initial slope (photosynthesis gained per unit of light when light is low), and \(R_d\) is the respiration rate (the leaf’s net CO₂ loss in the dark). In R: A ~ Amax * (1 - exp(-alpha * PAR / Amax)) - Rd.

Q8. Plot A against PAR for the sun leaves. From your plot, write down starting values for Rd (minus the value of A at PAR = 0), Amax (the level where the curve flattens, plus Rd) and alpha (the rise in A over the first 200 µmol of PAR, divided by 200). Fit the model to the sun leaves with nls() and report the three estimates.

sun <- subset(leaf, leaf_type == "sun"); nls(A ~ ..., data = sun, start = list(Amax = ..., alpha = ..., Rd = ...)). tapply(sun$A, sun$PAR, mean) helps you read the plot. Script Part F, sections 30 and 33; slide nls() needs starting values.

Q9. Report Amax and alpha for the sun leaves with 95% confidence intervals and write one results sentence, with units.

confint(model). Script Part F, section 31; slide A confidence interval for Nc.

Q10. Fit the same model to the shade leaves. Plot both fitted curves with their 95% bands on the data. Do sun and shade leaves differ in Amax? Answer with the confidence intervals, and suggest one biological reason.

nls_ci(model, data.frame(PAR = seq(0, 2000, 10))) from the code page gives fit, lwr and upr for the band. Script Part F, sections 31 and 32; slides A confidence interval for Nc and The quadratic overshoots.